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In  this  paper  results  are  presented  of  a  numerical  study  performed  for  four  different  burner  configura¬ 
tions  in  a  furnace  equipped  with  three  pairs  of  flameless  combustion  burners  firing  Dutch  natural  gas. 
The  simulations  have  been  validated  against  previously  published  results  of  an  experimental  study  [1], 
The  commercial  Computational  Fluid  Dynamics  (CFD)  code  Fluent  6.3  was  used  for  the  calculations. 
Using  the  Eddy  Dissipation  Concept  (EDC)  model  for  turbulence-chemistry  interaction  in  combination 
with  the  realizable  k—e  model  for  turbulence  and  a  skeletal  chemistry  mechanism,  the  main  furnace 
performance  was  consistently  reproduced  for  all  the  investigated  burner  configurations.  Moreover,  it  was 
found  that  due  to  relatively  low  Reynolds  numbers  in  the  cooling  air  flow  in  the  annulus  of  the  cooling 
tubes,  predictions  of  the  heat  extraction  rates  of  these  cooling  tubes  were  improved  by  treating  the  flow 
in  the  cooling  tubes  as  laminar.  Furthermore,  the  applied  error  tolerance  of  the  1SAT  procedure  was 
insufficient  for  accurate  species  concentration  predictions,  however,  based  on  analysis  of  the  main 
species  concentrations  in  the  flue  gas,  this  inaccuracy  did  not  influence  the  overall  predictions. 

The  most  important  experimental  results  have  been  investigated  using  the  CFD  simulations.  Firstly, 
a  longer  path  length  from  the  firing  burners  to  the  stack,  compared  to  the  path  length  to  the  regenerating 
burners,  explained  the  lower  CO  emissions  in  the  flue  gas  in  the  stack.  Secondly,  it  was  found  that 
a  recirculation  zone  between  the  upper  firing  burners  and  the  stack  in  configurations  C4  and  C5  resulted 
in  a  smaller  fraction  of  the  flue  gases  leaving  the  furnace  via  the  stack  compared  to  the  other  configu¬ 
rations.  Thus,  a  larger  fraction  left  the  furnace  via  the  regenerating  burners  and  this  resulted  in  higher 
preheat  temperatures  of  the  combustion  air.  Furthermore,  more  pronounced  recirculation  zones  in 
configurations  C3  and  C4  led  to  higher  temperature  uniformities  in  the  furnace.  Finally,  it  was  confirmed 
that  the  jets  of  the  burners  in  configurations  Cl  and  C3  showed  similar  merging  behavior,  leading  to 
similar  NO  emissions,  as  observed  in  the  experiments. 

©  2011  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

The  industrial  wish  for  higher  energy  efficiency  of  large-scale 
furnaces,  and  its  associated  fuel  savings,  has  demanded  for  new 
combustion  technologies,  combining  heat  recirculation  from  the 
flue  gas  with  low  pollutant  emissions.  Flameless  combustion  (also 
known  as  flameless  oxidation  (FLOX)  [2],  HiTAC  [3]  or  MILD 
combustion  [4])  is  such  a  novel  combustion  technique.  In  this 
technique  the  combustion  air  can  be  highly  preheated,  without 
increasing  the  pollutant  emissions,  in  particular  NO*.  The  air  and 
fuel  are  injected  at  high  velocity  and  spatially  separated.  Due  to  this 
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high  momentum  injection,  large  quantities  of  hot  flue  gases  are 
entrained  into  the  jets  before  they  mix  with  each  other,  lowering 
the  oxygen  availability  in  the  reaction  zone,  thus,  also  lowering  the 
local  reaction  rates  and  the  peak  temperatures.  These  low  peak 
temperatures  reduce  the  thermal  NO  emissions  [5], 

Since  the  introduction  of  flameless  combustion  in  the  early 
nineties  of  the  last  century,  many  universities  and  research 
departments  of  industry  have  made  efforts  in  experimentally 
investigating  this  new  technology.  These  studies  have  been  per¬ 
formed  on  many  different  scales,  from  single  open  flames  in  jet-in- 
hot-coflow  setups  [6,7]  to  (semi-)industrial  scale  test  furnaces 
equipped  with  multiple  burners  [8—11  ].  A  focus  is  here  on  furnaces 
equipped  with  multiple  flameless  combustion  burners. 

Many  of  these  experimental  studies  have  been  complemented 
with  Computational  Fluid  Dynamics  (CFD)  simulations.  In  the 
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Fig.  1.  Furnace  sketch,  representing  burner  configuration  C5  firing  in  parallel  mode.  The  boxed  numbers  1  and  2  indicate  the  two  sample  positions  for  the  flue  gas.  Sampling  point  2 
is  after  the  regenerators.  The  vertical  symmetry  plane  is  indicated  by  the  (green)  shaded  plane.  All  dimensions  are  in  mm.  The  two  inserts  show  enlarged  front  views  of  the  mesh 
around  a  burner  (left  hand  side)  and  a  cooling  tube  (right  hand  side).  The  total  mesh  contains  approximately  1.5  million  hexahedral  cells.  (For  interpretation  of  the  references  to 
colour  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


following  the  quality  of  these  simulations  of  furnaces  with  multiple 
flameless  combustion  burners  and  the  choice  of  physical  models  in 
these  simulations  are  discussed. 

The  200  l<W_th  semi-industrial  furnace  at  the  Royal  Institute  of 
Technology  (KTH)  has  been  simulated  extensively  using  the  STAR- 
CD  CFD  package  [12].  Besides  the  standard  k—e  turbulence  model, 
several  combustion  models  have  been  investigated,  of  which  the 
combination  of  the  Eddy  Dissipation  model  and  Finite  Rate  chem¬ 
istry  (ED/FR)  turned  out  to  give  better  results.  A  two-step  chemical 
mechanism  for  the  combustion  of  the  fuel  (LPG,  i.e„  propane),  with 
CO  as  an  intermediate,  was  used.  Thermal  radiation  was  calculated 
with  the  Discrete  Transfer  method  and  the  absorption  coefficient  of 
the  gases  using  the  Weighted  Sum  of  Gray  Gases  (WSGG)  model. 
The  cooling  tubes  were  not  incorporated  in  the  simulation,  but 
a  temperature,  based  on  measured  temperatures,  was  set  as 
a  boundary  condition.  For  the  validation  of  the  simulations 
measured  heat  fluxes,  wall  temperatures  and  in-furnace  species 
concentrations  were  used  [13].  The  validation  results  were 
reasonable,  however,  the  in-furnace  temperatures,  which  this  set  of 
models  is  known  to  overpredict,  were  not  compared.  The  three 
experimentally  investigated  firing  modes,  i.e.,  which  combination 
of  burners  form  a  burner  pair,  were  also  investigated  numerically. 
The  main  differences  observed  in  the  experiments  were  repro¬ 
duced  by  the  simulations.  Initially,  the  analysis  of  the  results  was 
focused  on  the  fundamental  properties  of  flameless  combustion.  In 
a  later  study  (firing  natural  gas)  also  the  observed  and  reproduced 
differences  between  the  firing  modes  were  investigated  numeri¬ 
cally  in  more  detail  [14]. 


In  2004  Hekkens  et  al.  simulated  the  experiments  performed 
with  the  1000  kW_th  IFRF  furnace  firing  natural  gas  using  the 
Fluent  CFD  package  [15,16].  Again,  the  standard  k—e  turbulence 
model  was  used.  Three  different  combustion  models  were  applied; 
two  Probability  Density  Function  (PDF)  methods,  assuming  chem¬ 
ical  equilibrium  and  laminar  flamelets,  and  the  ED/FR  model.  In  this 
last  model  a  two-step  chemistry  for  the  combustion  of  methane 
was  used.  Radiation  was  incorporated  using  the  Discrete  Ordinates 
method  and  the  absorption  coefficient  was  calculated  with  the 
domain-based  WSGG  model.  After  detailed  comparison  of  the 
numerical  and  experimental  results,  it  turned  out  that  the  ED/FR 
model  performed  best  regarding  the  species  concentration 
predictions  inside  the  furnace.  The  PDF  methods  were  unable  to 
correctly  predict  the  temperatures  inside  the  furnace,  which  can  be 
explained  by  the  fact  that  they  assume  (too)  fast  chemistry. 
However,  also  the  ED/FR  model  overpredicted  the  peak  tempera¬ 
tures  inside  the  furnace,  even  though  two  model  constants  have 
been  adjusted  for  flameless  combustion  purposes  using  previous 
IFRF  measurements. 

The  200  l<W_th  natural  gas  fired  flameless  combustion  furnace 
at  the  Faculte  Polytechnique  de  Mons  has  been  simulated  by  Lupant 
et  al.  using  the  Fluent  CFD  package  [10,17,18],  The  standard  k—e 
turbulence  model  is  applied,  with  a  model  constant  adjusted  for 
improved  prediction  of  the  spreading  rate  of  the  jets.  Combustion 
modeling  is  done  using  both  the  PDF  method  assuming  chemical 
equilibrium  and  the  combined  ED/FR  model.  Again,  as  in  the 
simulations  for  the  IFRF  furnace,  two  model  constants  in  the  ED/FR 
model  were  changed  in  order  to  improve  the  results.  A  one-step 
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Fig.  2.  Overview  of  simulated  burner  configurations.  The  two  rectangles  represent  the  two  side  walls  of  the  furnace  with  the  burner  flanges.  The  large  circles  represent  the  burner 
flanges,  if  filled  black  it  is  occupied  by  a  firing  burner,  if  filled  gray  it  is  occupied  by  a  regenerating  burner.  The  black  and  gray  circles  switch  after  a  period  of  30  s.  The  small  circles 
denote  the  position  of  the  cooling  tubes. 
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Table  1 

Composition  of  actual  and  simulation  fuel  for  Dutch  natural  gas  [26]. 


Component 

DNG 

Actual 

%-mass 

Simulation 

%-mass 

ch4 

69.97 

75.98 

c2h6 

4.63 

- 

c3h8 

0.90 

- 

C4H,o 

0.47 

- 

C5H,2 

0.16 

- 

c6h,4 

0.23 

- 

02 

0.02 

0.02 

n2 

21.52 

21.87 

co2 

2.10 

2.13 

Fig.  3.  Schematic  representation  of  a  single  cooling  tube.  All  dimensions  are  in  mm. 


this  standard  set  of  models,  a  wide  variety  of  models  for  turbulence, 
wall  functions  and  PDF  shapes  were  compared  [21],  The  simula¬ 
tions  are  validated  against  in-furnace  temperature  measurements 
[11],  The  agreement  is  reasonable  outside  the  combustion  zone. 
Inside  the  combustion  zone  the  temperatures  are  fairly  over¬ 
predicted,  as  was  observed  also  in  the  other  simulations  using  a  PDF 
method.  After  the  validation  of  the  simulation,  several  numerical 
investigations  were  performed  [19,20]. 

Generally,  it  can  be  concluded  that  combustion  models 
assuming  fast  chemistry  (e.g.  PDF  methods  assuming  chemical 
equilibrium)  are  not  applicable  for  the  simulation  of  flameless 
combustion.  Also,  a  careful  choice  of  the  chemical  mechanism  is 
important. 

At  Delft  University  of  Technology  (DUT)  a  300  l<W_th  furnace 
equipped  with  three  pairs  of  flameless  combustion  burners,  with 
the  unique  possibility  to  vary  the  positions  of  the  burners  in  the 
furnace,  has  been  investigated  experimentally  [1,22].  In  this  paper 
the  results  of  detailed  CFD  simulations  of  this  furnace  are  pre¬ 
sented.  The  objective  of  the  simulations  was  to  investigate  and 
explain  the  observed  trends  in  the  furnace  performance  for  the 
different  burner  configurations  in  the  furnace. 


2.  Experimental  setup 


chemistry  mechanism  has  been  used  for  the  combustion  of 
methane.  Radiative  heat  transfer  is  modeled  using  the  Discrete 
Ordinates  method  and  the  WSGG  model  was  used  for  the  calcula¬ 
tion  of  the  absorption  coefficient.  Validation  of  the  simulation  is 
performed  by  comparison  of  the  numerical  results  with  in-furnace 
radiative  heat  flux,  temperature  and  species  measurements. 
Reasonable  agreement  was  achieved  for  the  heat  fluxes,  however, 
all  simulations  overpredicted  the  temperature  inside  the  furnace. 
This  difference  was  smaller  in  the  simulation  where  the  model 
constants  in  the  ED/FR  model  were  adjusted.  Although  the 
maximum  concentrations  of  O2  were  predicted  rather  well,  the 
concentration  dropped  too  fast  after  its  maximum.  This  indicates 
that  the  combustion  reactions  are  too  fast,  which  also  explains  the 
prediction  of  the  temperatures  being  too  high.  Finally,  it  is  noted 
that  the  modification  of  the  model  constant  in  the  turbulence 
model  improved  the  prediction  of  the  delay  in  the  reaction. 

Several  regenerative  reheating  furnaces  at  the  NKK  steel 
company  firing  Low  Calorific  Value  gases  were  simulated  by  Ishii 
et  al.  using  the  Fluent  CFD  package  [11,19,20].  Again,  the  standard 
k—e  model  for  turbulence  is  used.  For  the  modeling  of  the 
combustion  the  PDF  method  assuming  chemical  equilibrium  is 
applied.  Radiation  is  incorporated  using  the  P-1  model  and  again 
the  WSGG  model  was  used  for  the  absorption  coefficient.  Besides 


The  furnace  is  equipped  with  three  pairs  of  regenerative 
flameless  burners.  The  burner  pairs  have  a  thermal  power  of 
100  kW_th  each,  thus  300  l<W_th  in  total.  The  furnace  has  inner 
dimensions  of  1500  x  1500  x  1850  mm  (length  x  width  x  height). 
The  insulation  consists  of  three  layers  of  ceramic  fire  bricks, 
together  300  mm  thick.  In  total  18  flanges  for  the  burners  are 
divided  over  two  opposite  sides  of  the  furnace  (nine  each).  In  this 
way,  it  is  possible  to  investigate  different  burner  configurations  in 
the  furnace.  The  heat  sink  is  made  of  eight  single  ended  cooling 
tubes,  four  placed  above  the  burners  and  four  beneath,  see  also 
Fig.  1  and  Section  3.2. 

The  burners  were  manufactured  by  Warmeprozesstechnik 
GmbH  and  are  of  the  REGEMAT  CD  200  type.  Each  burner  has  four 
combustion  air/flue  gas  nozzles  around  a  central  fuel  nozzle,  with 
diameters  of  20  mm  and  12  mm,  respectively.  The  burners  have 
ceramic  honeycomb  heat  exchangers  (regenerators)  incorporated. 
Eighty  percent  on  mass  basis  of  the  flue  gas  is  sucked  by  a  fan  via 
the  air  nozzles  over  these  honeycombs  (thus  heating  them)  for 
regeneration  of  the  heat,  while  the  remaining  twenty  percent  is 
leaving  the  furnace  via  the  central  stack  at  the  roof.  Three  burners 
are  firing  simultaneously,  while  the  other  three  burners  are 
regenerating.  After  a  preset  time  interval  of  10—30  s  they  switch 
and  the  firing  burners  start  regenerating,  and  vice  versa.  The  cycle 
time  (tcycie)  is  defined  as  one  periodic  cycle,  i.e.,  the  sum  of  the  firing 
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Fig.  4.  Predicted  velocity  magnitude  profiles  (m/s)  in  the  cooling  tube  at  three  axial  distances  z  (m)  versus  the  radial  distance  x  with  x  0  at  the  axis,  (a)  Realizable  k—e  turbulence 
model  and  (b)  laminar  zone. 
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Table  2 

Heat  extraction  rate  (Qtube)  and  temperature  of  the  cooling  air  at  the  outlet  (rcooiout) 
of  a  cooling  tube  from  configuration  Cl. 


Variable 

Experiment 

Turbulent  CFD 

Laminar  CFD 

Qtube  d<W) 

12.3 

18.9 

14.6 

^coolout  (°C) 

496 

702 

550 

Fig.  5.  Furnace  temperature  (r^m.  °C)  and  preheat  temperature  of  the  combustion  air 
(Fpreheat.  °C)  over  a  period  of  5  entire  cycles  from  the  C3  burner  configuration 
experiment. 

and  regenerating  period  for  one  burner.  The  experimental  data 
with  the  highest  cycle  time  (60  s)  were  taken  for  the  present  steady 
RANS  simulations,  since  these  experiments  are  considered  to  be 
closest  to  steady  state  operation.  The  data  with  an  excess  air  ratio 
(A)  of  1.25  were  used  for  the  simulations.  Dutch  natural  gas  (DNG) 
was  used  as  the  fuel,  which  has  a  net  (lower)  calorific  value  of 
31.669  MJ/m3.  The  furnace  operates  at  atmospheric  pressure.  Using 
an  ND1R  gas  analyzer,  the  NO  and  CO  concentrations  in  the  flue  gas 
from  the  stack  and  from  the  regenerators  (boxed  numbers  1  and  2 
in  Fig.  1,  respectively)  were  measured  on-line.  In  the  same  positions 
the  O2  concentration  was  determined  on-line  paramagnetically  for 
normalization  purposes.  More  details  on  the  experimental  setup 
and  results  can  be  found  in  previous  publications  [1,22]. 

Four  different  burner  configurations  were  investigated  numeri¬ 
cally.  For  all  of  these  configurations  the  firing  mode  was  parallel, 
indicating  that  all  the  firing  burners  are  on  one  side  of  the  furnace, 
while  all  the  regenerating  burners  are  on  the  other  side  of  the 
furnace.  In  Fig.  2  the  simulated  burner  configurations  are  presented. 

3.  Numerical  setup 

A  three-dimensional  mesh  of  the  furnace  was  generated  using 
Gambit  2.4,  see  also  Fig.  1.  The  mesh  contains  approximately  1.5 
million  hexahedral  cells.  Exploiting  the  symmetry  of  the  furnace 
with  respect  to  its  vertical  midplane,  to  reduce  the  computational 
time,  only  half  of  the  furnace  is  meshed.  The  cooling  tubes  are 
entirely  incorporated  in  the  simulations. 


The  simulations  were  performed  with  the  commercial  code 
Fluent  6.3.26  [23].  The  second  order  implicit  steady  solver  was  used 
with  Green-Gauss  cell  based  interpolation.  A  SIMPLE  algorithm  is 
used  for  pressure-velocity  coupling. 

The  calculations  were  performed  on  a  cluster  with  4  dual  core 
AMD  Opteron  8222  (3.0  GHz)  processors  with  8  GB  DRAM  per 
processor,  thus  32  GB  in  total.  The  simulations  were  calculated 
parallel  on  all  eight  nodes. 

3.2.  Furnace 

For  turbulence  closure  the  realizable  k—e  model  was  used.  This 
model  is  shown  to  give  improved  results,  compared  to  the  standard 
k—e  model,  for  the  prediction  of  round  jets  [24],  as  is  the  case  in 
these  burners. 

For  the  combustion  of  the  fuel  the  mechanism  of  Smooke  et  al. 
for  the  combustion  of  methane  was  applied  [25].  This  so-called 
skeletal  mechanism  consists  of  16  species  (resulting  in  16  trans¬ 
port  equations  to  be  solved)  and  incorporates  46  reactions.  It  is 
noted  here  that  a  mechanism  for  the  reaction  of  methane  only  is 
used.  This  implies  that  the  Dutch  natural  gas  (DNG)  is  replaced  by 
a  simulation  fuel.  This  simulation  fuel  contains,  as  a  fuel  compo¬ 
nent,  only  methane,  also  replacing  the  higher  hydrocarbons  present 
in  DNG,  in  such  way  that  the  heat  release  per  unit  mass  stays  equal. 
In  Table  1  the  composition  of  the  actual  DNG  and  its  simulation  fuel 
is  presented. 

The  turbulence— chemistry  interaction  is  taken  into  account  by 
using  the  Eddy  Dissipation  Concept  (EDC)  model  [27].  The  EDC 
model  has  been  successfully  applied  for  various  turbulent  and 
flameless  combustion  applications  [28,29].  Parente  et  al.  compared 
results  for  flameless  combustion  with  different  turbu¬ 
lence-chemistry  interaction  models  and  chemistry  mechanisms 
[30].  They  concluded  that  the  EDC  model  combined  with  a  skeletal 
or  full  chemistry  mechanism  performed  best.  Additionally,  Stefa- 
nidis  et  al.  also  compared  the  EDC  model  with  the  combined  ED/FR 
approach  [31  ].  They  report  that  the  main  disadvantage  of  the  ED/FR 
model  is  too  fast  reaction  rates  for  the  oxidization  of  the  fuel.  This 
disadvantage  is  especially  undesired  for  flameless  combustion 
applications,  since  in  flameless  combustion  the  reaction  rates  are 
relatively  low.  Finally,  Kim  et  al.  performed  a  study  on  the  use  of 
different  chemistry  mechanisms  for  the  simulation  of  flameless 
combustion.  The  EDC  model  was  used  in  all  the  test  cases,  since  it 
was  considered  by  the  authors  as  a  very  suitable  model  for  the 
simulation  of  flameless  combustion  [32]. 

The  in  situ  adaptive  tabulation  (ISAT)  procedure,  for  the  efficient 
computation  of  homogeneous  reactions  in  chemically  reacting 
flows,  has  been  applied  [33].  The  error  tolerance  of  this  procedure 
was  stepwise  decreased  until  a  value  of  1CT4.  For  one  case,  reported 
below,  the  error  tolerance  was  decreased  further  to  a  value  of  1CT5. 

Radiative  heat  transfer  is  accounted  for  using  the  Discrete 
Ordinates  method  and  the  domain-based  Weighted  Sum  of  Gray 
Gases  model  for  the  absorption  coefficient  of  the  gas  mixture.  The 


Table  3 

Comparison  of  main  variables  between  experiments  (EXP)  and  simulations  (CFD)  for  the  four  different  burner  configurations  (Cl,  C3,  C4  and  C5). 


Variable 

Cl 

C3 

C4 

C5 

EXP 

CFD 

EXP 

CFD 

EXP 

CFD 

EXP 

CFD 

Tfurn  (”C) 

1036 

1121 

1051 

1113 

1074 

1142 

1052 

1128 

Tstack  (°C) 

- 

1182 

- 

1192 

- 

1216 

- 

1195 

Qtool  (kW) 

99 

104 

97 

101 

98 

102 

90 

99 

Qloss  (kW) 

145 

110 

122 

108 

141 

110 

150 

110 

<&.  (kW) 

137 

- 

112 

- 

134 

- 

144 

- 

02  (%-vol,  dry) 

3.48 

3.22 

3.95 

4.76 

3.22 

3.47 

3.79 

3.43 

0|,h  (%-vol,  dry) 

3.66 

3.64 

5.58 

5.56 

4.16 

4.14 

4.03 

4.01 
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Fig.  6.  Main  predicted  results  for  configuration  C5  parallel  versus  ISAT  error  tolerance  (IET)  and  discretization  scheme;  ip3  =  1ET  of  10  3,  ip4  =  IET  of  10  4,  ip5  =  1ET  of  10  5, 
1  st  =  first  order  discretization  and  2nd  -  second  order  discretization,  (a)  Furnace  temperature  (Tfllrn,  °C)  and  flue  gas  temperature  in  the  stack  (rstack,  °C).  (b)  Heat  extraction  rate  of 
the  cooling  tubes  (Q,nn|,  kW)  and  other  heat  losses  (CLooi.  kW).  (c)  Oxygen  concentration  in  the  flue  gas  from  the  stack  (%-vol,  dry). 


mean  beam  length  was  based  on  the  dimensions  of  the  furnace 
since  the  flue  gas  composition  inside  the  furnace  is  reasonably 
uniform  and  with  this  approach  improved  predictions  of  the  radi¬ 
ation  fluxes  between  the  walls  and  cooling  tubes  are  expected.  The 
emissivity  of  steel  was  set  to  0.8  and  that  of  the  firebrick  insulation 
to  0.68. 

The  density  of  the  mixture  is  calculated  by  the  multi-component 
ideal  gas  law  using  the  low  Mach  number  approximation  p  =  pref. 
The  molecular  viscosity  and  thermal  conductivity  of  the  gas  mixture 
are  taken  constant  at  1.72  x  1CT5  kg/ms  and  0.0241  W/mK,  respec¬ 
tively.  The  mass  diffusivity  is  calculated  by  Fields  law  with  a  constant 
and  identical  mass  diffusion  coefficient  of  2.88  x  1CT5  m2/s  for  all 
species.  Finally,  the  heat  capacities  are  determined  per  species  with 
a  temperature  dependent  piece-wise  polynomial  provided  by 
Fluent,  and  for  the  mixture  with  the  mixing  law. 

The  heat  losses  through  the  walls  were  incorporated  by  calcu¬ 
lating  the  convective  and  radiative  heat  flux  to  the  environment. 
For  a  realistic  value  of  the  effective  thermal  conductivity  of  the 
walls,  the  predicted  wall  temperatures  were  compared  with  the 
experimental  values. 

In  previous  studies  this  set  of  models  proved  to  be  suitable  for  the 
simulation  of  fiameless  combustion  in  furnace  environments 
[34,14].  Additionally,  De  et  al.  performed  a  detailed  comparison 
study  of  these  and  comparable  models  for  experimental  results  of 
a  jet-in-hot-coflow  setup,  which  emulates  flameless  combustion 
[35],  Overall,  the  set  of  models  was  well  capable  of  predicting 
flameless  combustion.  Moreover,  it  was  shown  that  the  realizable 
k—t  performed  better  than  other  two-equation  turbulence  models. 
Also,  the  use  of  the  above-mentioned  constant  values  for  the 
viscosity  and  the  species  mass  diffusivities  showed  no  significant 
differences  with  simulations  with  a  temperature  dependent 
viscosity  or  using  multi-component  diffusion.  The  largest  difference, 
compared  to  the  experiments,  was  that  the  ignition  was  predicted  to 
occur  too  early.  Finally,  a  limitation  of  the  use  of  the  EDC  model  with 
a  skeletal  chemistry  mechanism  was  demonstrated  for  flows  with 
low  turbulence  Reynolds  numbers.  In  the  present  furnace  the  Rey¬ 
nolds  numbers  were  verified  to  be  higher  than  the  Reynolds 
numbers  where  these  limitations  become  apparent. 


H„0  46% 


o2 11% 

€ 

f  C02  43% 


0„  12% 


Ho0  45% 


CO<  1% 


CO„  42% 


CH,  <  1% 
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3.2.  Cooling  tubes 

In  preliminary  simulations  of  the  furnace  the  heat  extraction 
rate  of  the  cooling  tubes  was  fairly  overpredicted.  In  order  to 
investigate  this  difference  a  separate  simulation  was  performed  on 
a  single  cooling  tube.  The  mesh  of  this  single  cooling  tube  is  kept 
exactly  the  same  as  the  mesh  of  a  cooling  tube  in  the  mesh  of  the 
furnace.  As  a  boundary  condition  the  furnace  temperature  was  set 
on  the  outer  tube  in  this  separate  simulation. 

Each  cooling  tube  consists  of  two  concentric  annular  tubes.  Air 
was  used  as  the  cooling  medium  and  enters  the  inner  tube,  turns  at 
the  end  and  flows  back  through  the  annulus  between  the  inner  and 
outer  tubes.  This  setup  was  chosen  to  minimize  the  temperature 
gradients  along  the  length  of  the  outer  tube,  thus,  creating  an  as 
uniform  as  possible  heat  extraction  distribution.  The  inner  tube  has 
an  internal  diameter  of  100  mm  and  penetrates  1282  mm  into  the 
outer  tube,  which  in  turn  penetrates  1450  mm  into  the  furnace  and 
has  a  inner  diameter  of  150  mm.  In  Fig.  3  a  schematic  representa¬ 
tion  of  a  single  cooling  tube  is  presented. 

In  order  to  characterize  the  flow  regime  inside  the  cooling  tube 
the  Reynolds  number  was  calculated.  The  Reynolds  number  (Re)  is 
defined  in  Equation  (1). 


Re  = 


pudh 


(1) 


where  p  is  the  mean  density  (kg/m3),  u  the  mean  velocity  magni- 
tude  (m/s),  d/j  the  hydraulic  diameter  (m)  and  n  the  molecular 
viscosity  (kg/ms).  The  hydraulic  diameter  d/7  is  0.1  m  and  0.05  m  for 
the  inner  tube  and  the  annulus,  respectively.  The  Reynolds  number 


Fig.  7.  Main  species  mass  fractions  (-)  in  the  flue  gas  in  the  stack  with  burner 
configuration  C5  for  (a)  IET  =  10-4  and  (b)  IET  =  10-5. 


(X  0.1)  (X  0.01) 

Fig.  8.  Comparison  of  experimental  values  for  CO  concentration  (0.1  x  ppmv  @3%02, 
dry,  EXP)  and  predicted  values  from  the  simulations  (0.01  x  ppmv  @3%02,  dry,  CFD)  in 
the  flue  gas  from  the  stack  and  from  the  regenerators  for  all  four  burner  configurations 
firing  in  parallel  mode. 
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Fig.  9.  (a)  Sideview  of  the  furnace  showing  the  pathlines  from  the  nozzles  of  a  firing  burner  in  configuration  C5.  Pathlines  are  colored  by  the  velocity  magnitude  (m/s),  (b)  CO 
concentration  (ppmv  @3%02,  dry)  versus  distance  along  the  paths  from  the  firing  burner  to  the  opposing  regenerating  burner  (solid  line)  and  to  the  stack  (dashed  line).  Results  from 
CFD  simulations. 


of  the  flow  in  the  inner  tube  and  annulus  are  around  20,000  and 
10,000,  respectively.  These  values  indicate  that  the  flow  in  the  inner 
tube  is  turbulent,  however,  in  the  annulus  it  is  in  the  range  of 
Reynolds  numbers  that  indicate  the  transition  regime  between 
laminar  and  turbulent  flow. 

In  order  to  further  investigate  the  flow  regime  in  the  annulus,  the 
predicted  velocity  magnitude  profiles,  using  the  turbulence  model, 
in  the  cooling  tube  at  three  axial  positions  ( z )  are  presented  in 
Fig.  4(a).  Indeed,  the  velocity  profile  in  the  inner  tube  has  the  shape 
typical  for  turbulent  pipe  flow.  However,  the  velocity  profile  in  the 
annulus  tend  to  a  parabolic  form,  which  indicates  laminar  flow. 

In  a  subsequent  simulation  the  entire  cooling  tube  has  been 
treated  as  a  laminar  zone.  In  Table  2  the  numerical  results  for  the 
single  cooling  tube  (with  boundary  condition  values  taken  from  the 
Cl  configuration  experiment),  treating  the  flow  in  the  cooling  tube 
either  as  turbulent  or  as  laminar,  are  compared  with  the  experi¬ 
mental  values.  It  is  observed  that  the  predictions  for  both  the  total 
heat  extraction  rate  (Qmbe)  and  for  the  temperature  of  the  cooling  air 
at  the  outlet  (TCOoiout)  were  closer  to  the  experimental  values  when 
the  cooling  tubes  were  treated  as  laminar.  This  can  be  explained  by 
the  flow  in  the  annulus  to  rather  resemble  laminar  flow  behavior 
and  thus  creating  a  sort  of  thermal  resistance,  limiting  the  total  heat 
extraction  rate  of  the  cooling  tube.  In  Fig.  4(b)  the  velocity  profiles 
using  the  laminar  model  are  presented.  It  is  observed  that  the 
velocity  profiles  in  the  annulus  are  more  pronounced  parabolic. 

To  conclude,  in  the  furnace  simulations,  the  flow  in  the  cooling 
tubes  is  treated  separately  from  the  rest  of  the  furnace  and  as 
laminar  flow.  The  flow  within  the  furnace  itself  is  still  treated  as 
turbulent  flow. 


3.3.  Data  comparison 

For  comparison  of  the  transient  experimental  data  (due  to  the 
switching  of  the  burners)  and  the  steady  state  numerical  results, 
averaging  is  required.  In  Fig.  5  the  furnace  temperature,  charac¬ 
terized  by  a  thermocouple  in  the  furnace  side  wall,  and  the  preheat 
temperature  of  the  combustion  air  in  one  burner  are  presented  over 
a  period  of  5  entire  cycles.  These  data  are  from  the  C3  burner 
configuration  experiment. 

It  is  noted  that  the  transient  behavior  is  well  represented  in  the 
preheat  temperature,  while  the  furnace  temperature  is  not 
showing  any  correlation  with  the  cycle  time.  Actually,  all  the 
measurements  used  as  boundary  conditions  or  for  validation  of 
the  simulations  (wall  temperatures,  heat  extraction  rates  of  the 
cooling  tubes,  composition  of  the  flue  gas)  show  similar  steady 
signals. 

The  preheat  temperature  of  the  combustion  air  is  also  used  as 
a  boundary  condition  of  the  simulations,  since  the  regenerators  are 
not  incorporated  in  the  simulations.  Therefore,  the  data  of  the 
preheat  temperature  of  an  entire  firing  period  is  averaged  and  set  as 
the  boundary  condition. 

Furthermore,  in  order  to  have  a  similar  spatial  heat  distribution 
in  the  simulations  as  in  the  experiments,  two  separate  steady  state 
simulations  were  performed  for  each  configuration.  For  both  parts 
of  an  entire  cycle  a  steady  simulation  was  performed,  i.e.,  one 
simulation  with  one  set  of  burners  firing  and  one  simulation  with 
the  other  set  of  burners  firing  was  performed,  see  also  Fig.  2.  For  the 
final  comparison  of  the  experimental  and  numerical  results,  the 
results  of  these  two  simulations  were  averaged. 


a 

o 


Fig.  10.  (a)  Experimental  values  of  the  combustion  air  preheat  temperatures  (°C)  for  all  configurations,  parallel  firing  mode,  all  values  of  A  and  cycle  time  [1],  (b)  Position  of 
thermocouple  in  regenerators. 
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4.  Results  and  discussion 

4.3.  Validation 

In  Table  3  the  main  results  of  the  experiments  are  compared 
with  the  numerical  results  for  all  the  configurations.  It  is  noted  here 
that  all  the  four  configurations  are  simulated  with  the  exact  same 
set  of  physical  models  and  settings,  as  described  above. 

The  furnace  temperature  is  overpredicted  in  the  simulations 
with  around  70°  C,  however,  these  differences  are  considered  to  be 
acceptable. 

The  heat  extraction  by  the  cooling  tubes  (Qcooi)  is  predicted 
reasonably  accurately.  The  predicted  values  are  slightly  higher  for 
all  configurations,  but  the  values  are  close.  The  other  heat  losses 
(Qioss).  mainly  through  the  walls,  are  underpredicted  in  the  simu¬ 
lations.  This  is  partly  due  to  the  calculation  method  of  this  value  in 
the  experiments.  In  the  analysis  of  the  experimental  results,  the 
heat  flux  through  the  stack  is  calculated  assuming  that  the  flue  gas 
leaving  the  furnace  via  the  stack  is  at  a  temperature  equal  to  the 
furnace  temperature.  However,  comparing  the  furnace  tempera¬ 
ture  (Tfum)  and  the  flue  gas  temperature  in  the  stack  (Tstack)  in  the 
simulations,  the  latter  is  consistently  higher.  This  indicates  that  by 
taking  the  furnace  temperature  as  the  stack  flue  gas  temperature, 
the  experimental  heat  flux  through  the  stack  is  underestimated. 
The  values  for  the  experimental  heat  losses  calculated  with  the 
predicted  temperature  of  the  gas  in  the  stack  (Qj*oss)  show  that  this 
error  explains  only  part  of  the  difference. 

Finally,  the  oxygen  concentration  (O2)  in  the  flue  gas  is  consid¬ 
ered.  In  Table  3  the  measured  and  predicted  values  of  the  O2 
concentration  in  the  flue  gas  from  the  stack  are  presented.  Addi¬ 
tionally,  the  theoretical  value  (Ojf)  calculated  from  the  measured 
flows  of  comburants  for  complete  combustion,  both  using  the 
actual  fuel  (DNG)  and  using  the  simulation  fuel  (only  containing 
CH4  as  a  hydrocarbon,  see  above),  are  reported.  In  the  first  place, 
the  measured  values  are  compared  with  the  theoretical  values  for 
the  experiments.  All  the  measured  O2  values  are  lower  than  the 
theoretical  values.  These  differences  are,  besides  the  error  in  the 
measurement  of  the  oxygen  concentration  itself,  due  to  errors  in 
the  measurement  of  the  comburant  mass  flows.  Secondly, 
comparing  the  theoretical  values  for  the  actual  fuel  (DNG)  and  the 
simulation  fuel,  it  is  noted  that  they  are  very  similar.  Finally, 
comparing  the  predicted  values  in  the  simulations  with  the  theo¬ 
retical  values  for  the  simulation  fuel,  it  is  noted  that  the  predicted 
values  are  systematically  lower.  This  is  due  to  the  integration  of  the 
chemistry  in  the  simulations.  As  stated  before,  the  ISAT  procedure 
was  used  for  shorter  computational  times.  Gordon  et  al.  investi¬ 
gated  the  maximal  value  for  the  ISAT  error  tolerance  (IET)  for 
several  mechanisms  [36],  concluding  a  maximal  value  of  1CT6  for 
the  Smooke  mechanism  is  required.  Due  to  too  long  computational 
times,  this  value  was  not  supportable  for  the  current  mesh  in 
combination  with  the  current  cluster. 

To  verify  the  impact  of  the  error  due  to  this  deficiency,  one  case 
(configuration  C5  with  the  black  burners  firing,  see  Fig.  2)  was  run 
with  IET  1CT5.  In  Fig.  6  the  main  results,  as  presented  in  Table  3,  are 
plotted  for  the  simulations  with  the  different  values  of  the  IET.  It  is 
noted  that  a  large  difference  in  the  main  results  is  observed  between 
the  simulations  with  IET  1CT3  and  1CT4.  In  case  of  IET  1CT5  the 
results  do  improve,  however,  relatively  marginally.  It  is  concluded 
here,  that  these  simulations  would  improve  by  additional  calculation 
with  IET  1CT5,  however,  considering  the  time  investment  versus  the 
expected  improvement,  this  was  not  performed. 

Additionally,  in  Fig.  6,  the  results  of  the  simulation  with  IET  1CT4 
with  both  first  order  and  second  order  discretization  schemes  are 
compared.  It  is  noted  that  the  differences  in  main  results  between 
these  two  simulations  are  marginal. 


Fig.  11.  Comparison  of  the  maximum  value  of  t (-)  at  the  hot  side  of  the  regenerators 
in  the  experiments  (EXP)  and  the  predicted  r  (-)  of  the  flue  gas  entering  the  regen¬ 
erating  burners  (CFD). 

Finally,  in  Fig.  7  the  mass  fractions  of  the  main  species  in  the  flue 
gas  in  the  stack  are  compared  for  the  simulations  of  configuration 
C5  with  IET  1CT4  and  10~5.  It  is  noted  that  in  the  simulation  with  IET 
10  4  lower  O2  concentrations  and  higher  CO2  and  H2O  concentra¬ 
tions  are  predicted.  In  the  case  of  IET  1CT5  the  O2  content  is  slightly 
higher  and  the  CO2  and  H2O  contents  are  lower.  Thus,  by  decreasing 
the  IET  the  oxygen  concentration  is  tending  toward  the  experi¬ 
mental  and  theoretical  values.  For  both  cases,  the  concentrations  of 
CO  and  CH4  are  very  low;  they  are  in  the  order  of  10  and  1CT4  ppm, 
respectively.  It  can  be  concluded  that  in  the  simulation  with  the 
higher  value  of  IET  (1CT4)  the  combustion  is  complete,  i.e.,  there  is 
hardly  any  CO  and  CH4  left  in  the  flue  gas.  Therefore,  no  large  errors 
in  the  total  heat  distribution  in  the  furnace  due  to  integration  errors 
are  to  be  expected. 

Although  no  in-furnace  measurements  were  performed  for 
a  detailed  validation,  the  consistency  of  the  predictions  for  all  the 
four  investigated  configurations  confirmed  the  suitability  of  the  set 
of  physical  models  for  a  numerical  investigation  of  this  furnace. 

4.2.  Carbon  monoxide  emissions 

In  Fig.  8  the  concentrations  of  CO  in  the  flue  gas  from  the  stack 
and  from  the  regenerators  are  presented  for  the  experiments  and 
the  simulations.  For  these  two  flue  gas  sampling  locations,  see  the 
boxed  numbers  in  Fig.  1. 

All  the  predicted  values  of  the  CO  concentration  are  much  higher 
than  the  measured  values.  This  is  a  deficiency  of  the  used  chemistry 
mechanism  and  was  also  mentioned  as  such  by  the  developers  of  the 
mechanism  [25].  Additionally,  in  the  simulations  the  CO  concen¬ 
trations  in  the  flue  gas  from  the  regenerators  are  taken  before  the 
regenerators  (at  the  hot  side),  whereas  in  the  experiments  this  CO 
concentration  is  measured  after  the  regenerators  (at  the  cold  side). 
Thus,  part  of  the  excess  CO  observed  in  the  numerical  results  would 
have  been  converted  to  CO2  in  the  regenerators  if  simulated. 

However,  the  trend  in  the  difference  between  the  values  in  the 
flue  gas  from  the  stack  and  from  the  regenerators  is  reproduced 
correctly;  in  the  stack  flue  gas,  relatively  smaller  amounts  of  CO  are 
predicted,  whereas  in  the  regenerator  flue  gas  larger  amounts  of  CO 
are  predicted. 

The  difference  in  the  CO  concentrations  in  the  two  exits  is 
related  to  the  path  length  of  the  mean  flow  to  these  exits.  In 
Fig.  9(a)  the  pathlines  along  the  streamlines  of  the  mean  velocity 


Table  4 

Measured  mass  percentage  (%-mass)  of  flue  gas  leaving  the  furnace  via  the  stack. 


Cl 

C3 

C4 

C5 

Mass  percentage  stack  (%-mass) 

27.3 

31.2 

22.1 

22.1 
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(taken  from  the  CFD  simulations)  of  the  gases  flowing  from  the 
middle  firing  burner  in  configuration  C5,  to  the  opposing  regener¬ 
ating  burners  and  to  the  stack  are  presented.  The  burners  and  the 
stack  are  indicated  by  the  arrows.  It  is  observed  that  the  pathlines 
to  the  stack  are  longer  than  those  to  the  regenerating  burner.  Also, 
it  is  noted  that  part  of  the  gases  not  immediately  leaves  the  furnace, 
but  is  recirculated  to  the  root  of  the  comburants  jet. 

Based  on  the  observed  pathlines  from  CFD,  the  approximate 
path  lengths  of  the  mean  flow  were  calculated  exploiting  only 
straight  lines  and  right  angles.  The  recirculating  gases  were  not 
included  in  this  calculation  since  they  join  the  jet  again  at  its  root. 
Furthermore,  it  can  be  assumed  that  the  gases  that  immediately 
leave  the  furnace,  either  via  a  regenerating  burner  or  via  the  stack, 
account  for  the  largest  part  of  the  observed  CO  concentrations. 

The  path  lengths  between  all  firing  burners  and  all  regenerating 
burners  were  calculated  and  subsequently  averaged.  The  same  was 
done  for  the  path  lengths  between  the  firing  burners  and  the  stack. 
The  averaged  path  lengths  from  the  firing  burners  to  the  regener¬ 
ating  burners  and  to  the  stack  are  approximately  2.0  m  and  3.0  m, 
respectively.  For  the  firing  burner  in  Fig.  9(a)  these  path  lengths 


were  converted  to  residence  times  by  division  with  the  local 
velocity.  The  residence  time  of  the  gases  from  this  burner  to  the 
regenerating  burner  was  0.2  s.  From  the  firing  burner  to  the  stack 
the  residence  time  was  0.5  s.  Again,  these  values  include  the  direct 
paths  only,  excluding  recirculation. 

In  Fig.  9(b)  the  CO  concentration  along  the  considered  paths 
from  the  firing  burner  to  the  opposing  regenerating  burner  and  to 
the  stack  are  presented  from  the  simulations.  The  initial  1.5  m  is 
identical  for  both  paths,  since  in  both  cases  the  gases  cross  the 
furnace  width  first.  The  peak  in  CO  concentration,  which  can  be 
regarded  as  the  center  of  the  combustion  zone,  is  at  approximately 
1.2  m  from  the  burner  nozzles.  After  this  point,  the  CO  concentra¬ 
tion  decreases  steadily.  It  is  observed  that  the  CO  in  the  gases  to  the 
regenerating  burner  is  not  fully  converted  before  reaching  the  exit 
point.  On  the  other  hand,  in  the  gases  traveling  toward  the  stack, 
the  CO  concentration  further  decreases  (see  insert  with  enlarge¬ 
ment  in  Fig.  9(b)). 

Thus,  it  is  concluded  that  the  flue  gas  leaving  the  furnace  via  the 
stack  had  a  long  enough  residence  time  (due  to  its  longer  path 
length)  to  have  most  of  the  CO  converted  to  CO2. 


Cl  parallel 
Stack 


C5  parallel 
Stack 


Temperature 

(K) 

1 .706+03 
1 .636+03 
1  56e+03 
1 .496+03 
1 .426+03 
1 .356+03 
1 .286+03 
1.216+03 
I.He+03 
1.076+03 
1.006+03 
9.306+02 
8.60e+02 
7.906+02 
7  20e+02 
6  50e+02 
5.806+02 
5 10e+02 
4.409+02 
3  70e+0&7 
3.006+02 


C3  parallel 

Stack 


C4  parallel 


Stack 


Temperature 

(K) 

1.706+03 
1 .636+03 
1 .569+03 
1 .496+03 
1 ,42e+03 
1 .356+03 
1 .286+03 
1.218+03 
1.148+03 
1.07e+03 
1.006+03 
9.30e+02 
8.600+02 
7.90e+02 
7.20e+02 
6  500+02 
5.80e+02 
5.10e+02 
4.400+02 
3.700+02 
3.009+02 


Fig.  12.  Pathlines  from  all  grid  cells  of  the  firing  burner  inlets  of  all  configurations  colored  by  the  local  temperature  (K)  viewed  toward  the  vertical  midplane  of  the  furnace. 
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Fig.  13.  Cooling  tube  heat  extraction  rates  (kW)  comparison  for  the  experiments  (EXP) 
and  simulations  (CFD). 

4.3.  Combustion  air  preheat  temperature 

In  the  experiments  a  difference  in  the  preheat  temperature  of 
the  combustion  air  was  observed  for  the  different  burner  configu¬ 
rations.  In  Fig.  10  the  experimental  values  of  the  preheat  temper¬ 
ature  of  the  combustion  air  are  presented.  It  was  observed  that  the 
preheat  temperatures  for  the  C4  and  C5  configurations  are  higher 
than  those  for  configurations  Cl  and  C3. 

In  Fig.  11  the  maximum  value  of  the  normalized  temperature  (t) 
at  the  hot  side  of  the  regenerator  in  the  experiments  is  compared 
with  the  predicted  flue  gas  temperature  entering  the  regenerating 
burners.  The  normalized  temperature  t  is  defined  in  Equation  (2). 
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where  T  is  the  temperature  (K)  and  Tref  the  reference  temperature 
(273.15  K). 

The  numerical  values  are  higher  than  the  experimental  values. 
This  is  due  to  two  reasons.  In  the  first  place,  the  thermocouples 
inserted  in  the  regenerators  only  penetrate  until  the  end  of  the 
regenerators,  while  in  the  simulations  the  flue  gas  temperature  is 
taken  at  the  nozzle,  see  also  Fig.  10(b).  In  the  second  place,  the 
measured  temperatures  will  be  slightly  lower  than  the  actual  gas 
temperatures  due  to  radiative  heat  losses  of  the  thermocouples. 
The  reported  experimental  measurements  were  not  corrected  for 
these  radiative  losses. 

However,  in  both  the  experiments  and  the  simulations  the 
values  of  r  are  higher  for  configurations  C4  and  C5.  This  is  explained 
by  the  amount  of  flue  gas  leaving  the  furnace  via  the  stack.  In 
Table  4  the  experimental  values  of  the  mass  percentage  of  flue  gas 
leaving  the  furnace  via  the  stack  is  presented  for  the  four  config¬ 
urations.  It  is  noted  here  that  the  mass  percentage  of  flue  gas  to  the 
stack  is  set  as  a  boundary  condition  in  the  simulations. 

It  is  observed  that  for  both  Cl  and  C3  relatively  larger 
percentages  of  the  flue  gases  leave  the  furnace  via  the  stack.  This  is 


Fig.  14.  Temperature  uniformity  (-)  comparison  for  experiments  (EXP)  and  simula¬ 
tions  (CFD). 


Fig.  15.  Predicted  temperatures  (I<)  for  ail  thermocouples  in  burner  configurations  C4 
and  C5.  The  four  rectangles  represent  the  four  side  walls  of  the  furnace  with  the 
numbers  indicating  the  positions  of  the  thermocouples. 

due  to  the  fact  that  in  these  configurations  the  three  higher  burners 
are  positioned  in  the  upper  level,  whereas  in  C4  and  C5  they  are 
positioned  in  the  middle  level,  see  also  Fig.  2. 

In  Fig.  12  pathlines  along  the  streamlines  of  the  mean  velocity 
are  presented  from  all  grid  cells  of  the  firing  burner  inlets  of  all  the 
configurations,  colored  by  the  local  temperature.  It  is  observed  that 
in  configuration  C4  and  C5,  as  expected,  less  of  these  pathlines 
leave  the  furnace  via  the  stack  in  the  roof.  This  is  due  to  an  addi¬ 
tional  recirculation  zone  above  the  firing  burners  compared  to 
configuration  Cl  and  C3. 

It  is  concluded  that  in  configurations  Cl  and  C3,  due  to  the  firing 
burners  positioned  closer  to  the  stack  and  the  absence  of  a  recir¬ 
culation  zone  between  these  firing  burners  and  the  stack,  signifi¬ 
cantly  more  flue  gas  leaves  the  furnace  via  the  stack  and  relatively 
less  hot  flue  gas  is  available  for  regeneration,  leading  to  lower 
preheat  temperatures  of  the  combustion  air. 
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Fig.  16.  Cooling  tube  heat  extraction  rate  (kW)  versus  temperature  uniformity  (-)  for 
all  burner  configurations  in  parallel  firing  mode.  Black  symbols  are  the  experimental 
values,  gray  symbols  represent  the  numerical  values. 
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Fig.  17.  Development  of  the  velocity  in  the  z  -direction  (m/s),  temperature  (°C)  and  mixture  fraction  (?,  -)  for  the  lower  burners  in  configurations  Cl  and  C3  parallel.  The  data  are 
plotted  over  the  three  gray  lines  depicted  in  the  sketch  of  the  numerical  domain  of  configuration  C3.  The  (green)  shaded  plane  is  the  vertical  symmetry  plane.  (For  interpretation  of 
the  references  to  colour  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


4.4.  Cooling  tube  heat  extraction 

In  Fig.  13  the  experimental  and  numerical  values  of  the  cooling 
tube  heat  extraction  rate  of  the  four  configurations  are  compared. 
The  predicted  values  are  slightly  higher  than  the  experimental. 


However,  the  trend  between  the  different  cases  is  well  captured  by 
the  simulations. 

The  cooling  tube  heat  extraction  is  directly  related  to  the  total 
mass  flow  of  cooling  air  to  the  cooling  tubes.  In  the  total  mass  flow 
of  cooling  air  the  same  trend  is  observed  as  in  Fig.  13.  Since  the 
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measured  total  cooling  air  mass  flows  are  set  as  boundary  condi¬ 
tions  in  the  simulations,  the  trends  in  the  total  cooling  tube  heat 
extraction  rates  are  also  reproduced  in  the  simulations. 

4.5.  Temperature  uniformity 


In  Fig.  14  the  experimental  and  numerical  values  of  the  temper¬ 
ature  uniformity  of  the  four  configurations  are  compared.  The 
temperature  uniformity  ( Tu )  is  defined  as  in  the  following  equation. 


where  N  is  the  total  number  of  temperature  measurement  posi¬ 
tions,  T;  is  the  temperature  in  the  ith  position  and  T  is  the  mean  of 
all  the  temperature  measurements.  The  value  of  the  Tu  is  between 
0  and  1,  where  the  value  1  indicates  a  perfectly  uniform  furnace. 

The  trend  in  the  temperature  uniformities  agrees  well  between 
the  experiments  and  simulations.  Especially,  the  observation  that 
configurations  C3  and  C4  achieve  higher  temperature  uniformities, 
compared  to  configurations  Cl  and  C5,  is  reproduced.  However, 
there  is  an  opposing  trend  in  the  temperature  uniformities  of  C3 
versus  C4  and  Cl  versus  C5. 

In  Fig.  15  the  predicted  values  for  the  wall  temperatures  are 
presented  for  burner  configuration  C4  and  C5.  The  exact  positions 
of  the  thermocouples  are  denoted  in  the  drawings  of  the  four  side 
walls  of  the  furnace. 

It  can  be  observed,  comparing  the  thermocouples  in  the  highest 
level  (numbers  1,  3  and  5)  with  those  in  the  middle  or  lower  level 
(numbers  2,  4  and  6),  that  the  difference  in  the  temperature  is 
larger  for  configuration  C5,  which  results  in  a  lower  temperature 
uniformity.  Actually,  in  configuration  C4,  thermocouples  1,  3  and  5 
have  relatively  higher  values.  This  is  due  to  the  fact  that  a  stronger 
recirculation  zone  in  this  configuration  results  in  lower  tempera¬ 
tures  in  the  upper  zone  of  the  furnace,  see  also  Fig.  12.  This  also 
explains  the  lower  temperature  uniformity  in  configuration  Cl, 
where  the  least  pronounced  recirculation  zones  can  be  identified. 
For  the  values  of  the  thermocouples  in  the  opposite  wall  (numbers 
10  to  15)  the  same  behavior  can  be  observed. 

4.6.  Temperature  uniformity  versus  cooling  tube  heat  extraction 

In  the  experimental  results,  a  linear  correlation  was  observed 
between  the  cooling  tube  heat  extraction  and  temperature 
uniformity  [1].  In  Fig.  16  these  values  are  plotted  for  both  the 
experiments  and  the  simulations.  It  is  observed  that  the  decreasing 
trend  in  the  cooling  tube  heat  extraction  rate  with  increasing  value 
of  the  temperature  uniformity  is  not  reproduced.  Due  to  the  choice 
of  experiments  simulated  it  is  not  possible  to  either  confirm  or 
contradict  the  observations  of  the  experimental  results.  In  order  to 
be  able  to  do  this,  additional  simulations  of  cases  exploiting  more 
variation  in  the  values  of  the  cooling  tube  heat  extraction  and  the 
temperature  uniformity  are  required.  Due  to  time  limitations,  these 
additional  simulations  are  currently  not  performed. 

4.7.  Point  of  confluence 

In  a  previous  study  on  a  comparable  furnace  [14]  the  point  of 
confluence  of  the  comburant  jets  was  indicated  as  an  important 
factor  for  the  NO  formation.  In  this  study  the  NO  emissions  were 
similar  for  all  configurations  and  operating  conditions  [1].  To  verify 
whether  the  jets  show  similar  merging  behavior  in  this  furnace  the 
firing  burners  in  the  lowest  level  of  configurations  Cl  and  C3  are 
compared. 


In  Fig.  17  the  velocity  in  the  (axial)  z  -direction,  the  temperature 
and  the  mixture  fraction  are  plotted  on  lines  at  three  distances  from 
the  nozzles  of  these  burners.  The  distance  from  the  burner  nozzles  z 
is  normalized  by  the  diameter  of  the  air  nozzles  d  (0.02  m).  The  fuel 
nozzles  are  positioned  at  x  =  0.295  m  and  0.75  m  and  are  sur¬ 
rounded  by  two  high  momentum  air  nozzles.  The  mixture  fraction 
was  calculated  by  applying  Bilger’s  formula  [37], 

The  point  of  confluence  is  defined  as  the  point  where  the  low- 
momentum  jet  is  absorbed  by  the  high  momentum  jet(s)  [38]. 
The  further  the  point  of  confluence  from  the  nozzles,  the  more 
internal  flue  gas  recirculation  is  present  before  the  fuel  and  oxidizer 
jets  merge  and  lower  NO  emissions  are  expected. 

For  the  burner  at  x  =  0.295  m  no  clear  differences  can  be 
observed  between  configuration  C3  and  Cl,  both  seem  to  have 
a  point  of  confluence  at  z/d  =  6.  The  burner  atx  =  0.75  m  has  a  point 
of  confluence  at  z/d  =  8  in  both  configurations  (only  configuration 
C3  is  shown  in  Fig.  17,  but  the  burner  of  configuration  Cl  showed 
the  same  value).  Concluding,  it  is  confirmed  that  the  jets  of  the 
burners  in  configurations  Cl  and  C3  show  similar  merging 
behavior,  which  is  in  agreement  with  the  proposed  explanation  for 
the  similar  NO  concentrations  in  their  flue  gases.  Also,  no  distinct 
differences  in  the  maximum  temperatures  could  be  identified 
between  the  two  configurations.  Finally,  it  can  be  expected  that  the 
middle  burner  (at  x  =  0.75  m)  produces  less  NO  than  the  two  outer 
burners. 

5.  Conclusions 

A  numerical  study  has  been  performed  for  four  different  burner 
configurations  in  a  furnace  equipped  with  three  pairs  of  flameless 
combustion  burners.  A  carefully  selected  set  of  physical  models 
proved  to  be  adequate  to  reproduce  the  main  results  observed  in 
the  furnace  experiments.  The  simulations  have  been  validated 
against  the  main  results  of  the  furnace  experiments  and  by 
consistently  capturing  the  trends  in  these  results  for  the  different 
burner  configurations.  Moreover,  it  was  found  that  due  to  relatively 
low  Reynolds  numbers  in  the  cooling  air  flow  in  the  annulus  of  the 
cooling  tubes,  predictions  of  the  heat  extraction  rates  of  these 
cooling  tubes  were  improved  by  treating  the  flow  in  the  cooling 
tubes  as  laminar.  Furthermore,  the  applied  error  tolerance  of  the 
ISAT  procedure  was  shown  to  be  insufficient  for  accurate  species 
concentration  predictions,  however,  based  on  analysis  of  the  main 
species  concentrations  in  the  flue  gas,  this  inaccuracy  did  not 
influence  the  overall  predictions. 

The  most  important  experimental  results  have  been  investi¬ 
gated  using  the  simulations.  In  the  first  place,  a  longer  path  length 
from  the  firing  burners  to  the  stack,  compared  to  the  path  length  to 
the  regenerating  burners,  explained  the  lower  CO  emissions  in  the 
flue  gas  in  the  stack.  Secondly,  it  was  found  that  a  recirculation  zone 
between  the  upper  firing  burners  and  the  stack  in  configurations  C4 
and  C5  resulted  in  a  smaller  fraction  of  the  flue  gases  leaving  the 
furnace  via  the  stack  compared  to  the  other  configurations.  Thus, 
a  larger  fraction  left  the  furnace  via  the  regenerating  burners  and 
this  resulted  in  higher  preheat  temperatures  of  the  combustion  air. 
Furthermore,  the  total  cooling  tube  heat  extraction  was  directly 
related  to  the  total  amount  of  cooling  air.  In  the  fourth  place,  more 
pronounced  recirculation  zones  in  configurations  C3  and  C4  led  to 
higher  temperature  uniformities  in  the  furnace.  Finally,  it  was 
confirmed  that  the  jets  of  the  burners  in  configurations  Cl  and  C3 
showed  similar  merging  behavior,  leading  to  similar  NO  emissions, 
as  have  indeed  been  observed  in  the  experiments. 

It  can  be  concluded  that  by  the  numerical  study  more  insight  is 
attained  regarding  the  effects  of  burner  positioning  in  a  multi¬ 
burner  flameless  combustion  furnace.  In  the  future,  the  simula¬ 
tions  can  be  improved  by  further  decreasing  the  error  tolerance  of 
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the  ISAT  procedure.  Also,  in  order  to  verify  the  observation  in  the 
experiments  of  a  decreasing  cooling  tube  heat  extraction  with 
increasing  temperature  uniformity,  additional  simulations  are 
required. 
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